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ABSTRACT 

We have tested a previous analytical estimate of the dynamical friction timescale 
in Modified Newtonian Dynamics (MOND) with fully non-linear N-body simulations. 
The simulations confirm that the dynamical friction timescale is significantly shorter in 
MOND than in equivalent Newtonian systems, i.e. systems with the same phase-space 
distribution of baryons and additional dark matter. An apparent conflict between this 
result and the long timescales determined for bars to slow and mergers to be completed 
in previous N-body simulations of MOND systems is explained. The confirmation of 
the short dynamical-friction timescale in MOND underlines the challenge that the 
Fornax dwarf spheroidal poses to the viability of MOND. 
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1 INTRODUCTION 

Numerous observations indicate that galaxies and clusters of 
galaxies have gravitational fields that are stronger at large 
radii than the standard theory of gravity predicts from the 
conjecture that light is a fair tracer of mass. The standard 
interpretation of this phenomenon is that ~ | of the mass 
in the Universe is contributed by particles that are dark 
because they do not interact electromagnetically. The alter- 
native hypothesis, that the standard theory of gravity fails 
at low accel erations, was advance d by Milgrom (1983). Sub- 
sequently iBekenstein fc MilgromI (|l984l ) proposed a modi- 
fication of Poisson's equation that put the proposal on a 
quantitative basis. Recently interest in Modified Newtonian 
Dynamics (MOND) has increased with a growing awareness 
that relativistically covariant theories are possible that re- 
duce in the non-relativist ic limit to the Bekenstein-Milgrom 
theory (|Bekensteinll2004l '). 

One way to rule out MOND would be to detect the 
particles that are supposed to contribute most of the Uni- 
verse's mass, and several experiments are endeavouring to 
do this under the assumption that the particles have weak 
interactions. Another way to rule out MOND would be to 
show that it is inconsist ent with observatio ns of the growth 
of cosmic structure (e.g. lZlosnik et al]|2007h or the dynamics 
and evolution of galaxies. 

In standard gravity, dynamical friction (DF) is thought 
to play a major role in se veral astrophysical contexts (e.g. 
iBinnev fc Tremaind |2008| . §8.1). Ciotti & Binney (2004; 
hereafter CB04) showed that DF is a more potent phe- 
nomenon in MOND than in standard gravity, and argued 



that it would cause the globular clusters of dwarf spheroidal 
galaxies to spiral to the centres of their hosts within a dy- 
namical time. However, the results in CB04 were obtained 
under rather restrictive assumptions and by a non-standard 
argument that involves the fluctuation-dissipation theorem 
(see Bekenstein fc M aoz 1992). In particular, CB04 gave 
results for plane-parallel systems that are in the "deep- 
MOND" regime. Given the prediction of CB04 that DF is 
much more powerful in MOND than in Newtonian gravity, 
it is important to test their results with a different approach 
and in less specialized contexts. 

In this paper we use N-body experiments to explore DF 
in regimes in which departures from Newtonian gravity are 
of varying magnitude. The paper is organized as follows. In 
Section 2 we outline our methodology explain our choice of 
problem. Section 3 gives details of the code and the initial 
conditions of the simulations. Section 4 presents the results 
and Section 5 relates them to other recent work. The conclu- 
sions and their astronomical implications are summarized in 
Section 6. 



2 METHODOLOGY 

Two fundamental features of MOND preclude a direct ex- 
tension of the standard Chandrasckhar-Spitzer derivation of 
DF (Chandrasekhar 1942; S pitzer 1987. ): in MOND (i) forces 
are not additive and (ii) all two-body orbits are bound, so 
individual encounters do not have a finite duration. On the 
other hand, there is no fundamental obstacle to direct simu- 
lations of DF in MOND: the drag experienced by a massive 
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object as it moves through a swarm of less massive objects 
can be measured in an N-body simulation. 

In the classical approach to DF one imagines the back- 
ground swarm to be homogeneous and perpetrates the Jeans 
swindle to neglect the mean gravitational field of the swarm. 
Even in Newtonian gravity this step is of questionable va- 
lidity, and it is inadmissible in MOND because (i) the long- 
range nature of two-particle interactions in MOND implies 
that distant encounters are liable to dominate (in the New- 
tonian case they nearly do) and (ii) the non-linearity of 
the Bekenstein-Milgrom field equation implies that the im- 
pact of an encounter depends on the nature of the mean 
field. Therefore, any attempt to extend to MOND the clas- 
sical approach to DF leads to the consideration of the drag 
that a massive body experiences as it moves through a self- 
gravitating system of finite size. 

In order to compare the action of DF in MOND and 
Newtonian systems as cleanly as possible, we compare DF 
in MOND with DF in the equivalent Newtonian system 
(ENS); that is the Newtonian system in which the visible 
matter has exactly the same phase-space distribution as 
in the MOND system (M ilgrom 2001; Nipoti ot al. 2007tj; 
iNipoti. Londrillo fc Ciottil l2007d '). In the ENS the visible 
matter is enveloped in a dark-matter halo that contributes 
to the DF experienced by a massive body. 

There are two astrophysically important prob- 
lems for which DF is crucial: (i) the inward spi- 
ralling of a small system that has fallen i nto a 
larger one (e.g. ^■emaine. Ostriker fc SpitzeJ 19751: 



Table 1. Parameters of the simulations. 



Bontekoe fc van Albadal 19871 : iHernguist &c Weinberel 
19891 : I Arena fc BertinI 120071 ') . and (u) the slowing of a 



massive b a r by t he material in which it i s emb e dded (e.g. 
Weinberg! Il985l : iDebattista fc SellwoodI 1 19981 : ISellwoodI 
20061 . and references th erein). Our MOND-enabled N - 
body code (N-MODY: ICiotti. Londrillo fc Nipotil [ioo 



iNipoti. Londrillo fc Ciottil2007al : iLondriUo fc Nipotil 120081 'i 
is better adapted for problem (ii) because it uses a polar 
grid and therefore has higher resolution at small radii than 
at large. So we study the effect that MOND has on the rate 
at which the rotation of a rigid bar is slowed by DF. 



3 THE NUMERICAL SIMULATIONS 



With iBekenstein fc MilgromI l|l984) we assume that Pois- 
son's equation V^i^'^ = A-kGp should be replaced by the 
non-relativistic field equation 



IV0II 
ao 



AixGp, 



(1) 



where ||...|| is the standard Euclidean norm, and 4> is the 
gravitational potential for MOND; the gravitational accel- 
eration is g = — V(/> just as the Newtonian acceleration 
is g = — V</> . For a system of finite mass, V0 — > as 
||x|| ^ oo. The function ii{y) is constrained by the theory 
only to the extent that it must run smoothly from fj,{y) ~ y 
at ?/ ^ 1 (the so-called deep-MOND regime) to p(y) ~ 1 at 
y ^ 1 (the Newtonian regime), with the transition taking 
place at J/ « 1, i.e., when || Vi^iH is of order the characteristic 
acceleration ao — 1.2 x 10~^°ms~^ . In th e present work we 
adopt p.{y) = y/^Jl + y"^ (|Milgromlll983h . 



Name 


Gravity 


7 


MnM /M-t 

^ J_J iVl / ^ 




/ ^* 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


{7) 


MOO 


MOND 








1 


0.61 


8 


EGO 


Newton 





6.6 


1 


0.61 


16 


MOl 


MOND 








0.1 


1.06 


8 


EOl 


Newton 





23.1 


0.1 


1.06 


16 


M02 


MOND 








0.01 


1.88 


8 


E02 


Newton 





75.3 


0.01 


1.88 


16 


MIO 


MOND 


1 





1 


0.74 


8 


ElO 


Newton 


1 


7.0 


1 


0.74 


16 


Mil 


MOND 


1 





0.1 


1.27 


8 


Ell 


Newton 


1 


24.4 


0.1 


1.27 


16 


M12 


MOND 


1 





0.01 


2.26 


8 


E12 


Newton 


1 


79.2 


0.01 


2.26 


16 



(1): name of the simulation. (2): gravity law. (3) inner 
logarithmic slope of the stellar density distribution 
[equation JJJl]. (4): dark matter to baryonic mass ratio. (5): 
internal acceleration ratio k = GM,/(aor^). (6): initial angular 
frequency of the bar, in units of 0« = l/t*. (7): total number of 
particles in units of 10^. 

From Poisson's equation and equation ([l| it follows that 
the MOND and Newtonian gravitational accelerations are 
related by fi{g/ao)g = g''^ -I- S, where g = |jg||, and S is 
a solenoidal field dependent on the specific p considered. 
Since in general S 7^ 0, standard Poisson solvers cannot be 
used to develop MOND N-body codes; equation ([1} must be 
solved at each time step (Brada fc Milgrom 1999, Nipoti et 
al. 2007a). 



3.1 The code 

N-MODY is a parallel, three-dimensional particle-mesh code 
that can be used to run eith er MOND or Newtonian sim- 
ulations (Nipoti et al. 2007a, iLondrillo fc Nipotil lioosi ') . In 
the present study the spherical grid has 128 radial nodes, 
64 nodes in colatitude 9 and 128 nodes in azimuth (j), 
and the total number of particles is in the range A'^part ~ 
8 X lO'' - 1.6 X 10^ We verified with convergence exper- 
iments that these numbers of particles and grid points are 
sufficient to exclude that our results are significantly affected 

by discreteness effects. 

In previous papers (jCiotti. Nipoti fc LondriUol |2007| . 
Nipoti et al. 2007ac) we have used N-MODY to demonstrate 
significant differences in the operation of violent relaxation 
in systems in which MOND is important and in their New- 
to nian equivalents. We refe r readers to these papers and 
to ILondrillo fc Nipotil (|2008l ) for details of the code and its 
tests. 



3.2 Initial conditions 

The baryonic component of the initial conditi ons of the sim- 
ulations is described by a spherical 7-model l|Dehnenlll993l : 
[Trcmainc ct al. 1994) with density distribution 



p*(r) 



47r rT(r, + r)-*-7 ' 



(2) 
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where M« is the total stellar mass, r, is the scale radius, 7 is 
the inner l ogarithmic slope , and we consider the cases 7 = 
and 7 = 1 ()Hernauistl ll990h. For a given model, the MOND 
potential is easily calculated, and to each MOND model with 
potential corresponds an ENS with ^J^t — 4>, thus having 
a total density ptot{r) = V^0(r)/47rG'. In principle, such a 
distribution would have infinite mass, so we truncate it at 
r ~ lOr, . 

The particles of the stellar component are distributed in 
phase-space with the standard rejection technique, restrict- 
ing for simplicity to the fully isotropic case. The MOND 
distribution function i s obtained numerically w ith an Ed- 
dington inversion fe.g. iBinnev fc Tremainj|2008l . §4.3.1) 



fM{E) = 



1 



d 



dp, d<j) 



(3) 



where the upper integration limit reflects the far-fleld loga- 
rithmic behaviour of the MOND potential. 

In the equivalent Newtonian models, at variance with 
Nipoti et al. (2007c), we do not distinguish between dark 
and stellar particles, but we consider a single component 
distributed according to the numerical isotropic distribution 
function 



fN{E) 



d 



x/87r2 dE 



dpt. 
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(4) 



Note that the upper limit of integration is due to the finite 
mass of the system (a consequence of the density truncation 
at lOr.). 

The physical scales of the problem are introduced as fol- 
lows. First we identify each MOND initial condition by fixing 
a value for the dimensionless internal acceleration parame- 
ter K = GMf/{aQr1), so that Af» and r* are not indepen- 
dent quantities: in physical units, r, ~ 3.4k~^^^M^'^jq kpc. 
Mi. Mq. The time and velocity units are 



where M,,io 

t* = y/rf/GMl ~ 29.7k-^/''M,^_(o Myr, and = r^jt, ~ 
112/s:^/*M^^_/okms"^ (Nipoti et al. 2007a). The simulations 
are evolved up to f = SOTo (where To is the initial ro- 
tation period of the bar, see Section I3.3|l . with timestep 
= O.OOSTo. 

We verified that both the MOND system and its ENS 
were in equilibrium by ruiming them without the bar for 
several dynamical times. The dark and luminous matter 
distributions of the considered ENSs are shown in Fig. [TJ 
plotting, as functions of radius, their dark-matter density 
Pdm = Ptot — p* (top panel) and their total (dark-matter 
plus stars) mass-to-light ratio MjL = T -^ptot/ p*, where 
T» is the stellar mass-to-light ratio (bottom panel). All 
the considered ENSs are dark-matter dominated, because 
the corresponding MOND cases have internal accelerations 
<ao. Focusing on the central regions of the ENSs, we note 
that, for fixed k, models with cored stellar profile (7 = 0, 
dashed curves in Fig. [T} have higher M/L and shallower 
inner dark-matter profile than models with cuspy stellar 
profile (7 = 1, solid curves in Fig. In particular, we 
have 7dm ~ 1 when 7 = 1 and 7dm ~ 1/2 when 7 = 0, 
where 7dm = — limr-,o(rflnpdm/rflnr) is the inner logarith- 
mic slope of the dark-matter density distribution. 
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Figure 1. Dark-matter density pdm = ptot — p» (top panel) and 
total (dark-matter plus stars) mass-to-light ratio M/L (bottom 
panel) as functions of radius for the equivalent Newtonian models. 
M/L is in units of the stellar mass-to-light ratio T* . The first digit 
after "E" is the value of 7 while the second digit is — log]^Q(K) 
(Table 1). Dashed and solid curves refer to models with 7 = 
and 7 = 1, respectively. 



3.3 Bar equation of motion and the warm-up of 
the gravitational field 

The density distribution of the rigid bar is a prolate ellipsoid 
of density distribution 

Mb 



Pb 



where in the Cartesian system co-rotating with the bar 



,2 ^ y" 



g> 1. 



(5) 



(6) 



In all the cases here presented we use Mb = 0.05M*, 
b = 0.2r» and q = 5: with this choice, the length of the 
bar's semi-major axis is r*. We adopted the distribution ([Sjl 
because the Gaussian dependence on m produces a bar well 
defined spatially, but at the same time it avoids a density 
discontinuity that could produce spurious behaviour in the 
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Figure 2. Time evolution of the bar angular frequency in units of = l/t» in MOND (solid curves) and equivalent Newtonian models 
with dark matter (dotted curves), for six pairs of models with different combinations of the parameters k and 7. Cuspy models are on 
the right. Smaller values of k correspond to deeper MOND regimes. Time is normalized to the initial rotation period of the bar Tq. 



potential solver at the bar's edge. In all the simulations the 
bar rotates around the z' axis. 

The Lagrangian formulation can be used to derive the 
bar's equations of motion in an inertial Cartesian frame 
{x,y,z) such that z = z': with iph the angle in the {x,y) 
plane between the x and x' axes, we have 



J 33- 



p,{m)[y-^x-]dV, 



(7) 



where 

2 , . x2 2 (w cos G9b — a; sin o9b)'^ 
m = (a;cos((3b + 2/smiy3b) + z H (8) 

and from equation ((5} 

^ 4^g(l + g^) r i , Mbb^(l + g^) 

/33 = J m pb(m)dm = ^- -. (9) 

In order to start the numerical simulations smoothly, we 
take the density of the bar to be the product of equation ((5]) 
and the function 



a{t) = min ( 1 



(10) 



where the growth time is tg = 3To with To = 27r/S7b,o and 
S7b,o is the initial angular velocity of the bar's frame. 

As a consequence, when the bar mass is very small, the 
underlying density distribution is only slightly affected by 
its gravitational field, and the integral at the r.h.s of equa- 
tion (O almost vanishes. Of course, in the MOND simula- 
tions, on the r.h.s. of equation ([T} we have p — p, + pb- In all 
the presented cases we use Qh,o = Uc(r'*)/r*, where Uc(r'«) is 
the circular velocity at r* (i.e. at the edge of the bar). 

For all the simulations, the inner logarithmic slope 7, 
the internal acceleration ratio k, the total dark to luminous 
mass ratio Mum/M^, the initial angular velocity of the bar 
Qh,o and the total number of particles A^'part are given in 
Table 1. 
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4 RESULTS 

CB04 found that the DF timescales in a (deep) MOND sys- 
tem (tMc), and in the equivalent Newtonian system (tf^.^^) 
are related by 



where 1 + TZ = g/g^ with g the modulus of the MOND field 
and the modulus of the Newtonian field generated by the 
baryonic distribution. We recall again that the treatment 
of CB04 was carried out for a plane-parallel distribution 
of field particles. In the present case the situation is more 
complicate, as TZ is not constant in the system. Thus, for 
each pair of simulations with the same value of k and 7 (and 
therefore with the same gravitational field) we assume as 
fiducial value 1-1-7?. = gi(r«)/(;J^(r*), i.e. we evaluate TZ at the 
edge of the bar. Models with lower values of the parameter k 
are in deeper-MOND regime, so TZ increases for decreasing k. 
In particular, we considered the cases n — 0.01, k — 0.1 and 
K = 1, so all the models presented have internal accelerations 

For each simulation. Fig. [5] plots against time the an- 
gular frequency of the bar fib = V^b; models with constant- 
density cores are on the left and those with cuspy cores are 
on the right. Systems with smaller k and therefore more 
MOND/dark- matter dominated dynamics are lower down. 
In each panel the full curve shows the MOND system, and 
the dotted line its Newtonian equivalent. In every case the 
bar slows down more in the MOND system than in its ENS. 

The instantaneous DF timescale is | f2b /fib | , but in prac- 
tice this quantity fluctuates strongly, so is is not a useful 
measure of the DF timescale. A more robust measure is 
the time t^ic for fib to reach 70% of its initial value, so 
fib(ifric) = 0.7fib,o- FigureOis a plot against 1 -|- 7? of the 
MOND value of t^c divided by tfHc of the ENS. This ratio 
decreases with k as predicted by CB04, and in the deep- 
MOND regime the rate of decrease with increasing 1 + TZ 
is consistent with that predicted by CB04 (shown by the 
full line). The vertical offset between the full line and the 
points could partly reflect the arbitrariness of our defini- 
tion of tfric and our assignment of values of TZ to simula- 
tions. However, we verified that for any sensible definition 
of TZ, the experimental ratio tf^^^/tf^.^^ is always closer to 
unity than CB04 predicted. Presumably this finding arises 
because we are considering spherical systems, while CB04 
assumed plane-parallel symmetry. 



5 RELATION TO OTHER WORK 

iTiret fc Combej l|2007ah used N-body simulations similar, 
in some respect, to those reported here to compare the dy- 
namics of galactic bars in MOND and in equivalent New- 
tonian systems. They found that DF had a much smaller 
effect on bars in MOND than in the Newtonian cases. This 
finding is not in conflict with our opposed finding for this 
reason: whereas our bars contain only five percent of the 
galajcy's baryonic mass, in the Tiret & Combes simulations 
the bars contain the majority of the baryonic mass. Conse- 
quently, there is significant background mass for the Tiret 

6 Combes bars to interact with through DF only when dark 
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log(l+i?) 

Figure 3. Ratio of MOND and Newtonian equivalent DF 
timescale tf^ic as a function of the parameter Ti for pairs of mod- 
els with 7 = 1 (solid symbols) and 7 = (empty symbols). We 
measure tf^ic in the simulations as the time for f2b to reach 70% of 
its initial value. The solid line is the ratio as computed by CB04 
in the deep-MOND limit. 

matter is present. By contrast our bars interact with a back- 
ground that contains 95 percent of the baryonic mass, so DF 
is effective. 

Similar reason ing applies to the find ings of Nipoti 
et al. (2007c) and iTnet fc Combed (|2007bh that merging 
timescales are much longer in MOND than in Newtonian 
gravity with dark matter: in the Newtonian case there is 
an abundance of dark matter to absorb energy and angular 
momentum from the stellar systems, while in MOND the 
energy and angular momentum has to be absorbed by the 
outer parts of the stellar systems themselves. Consequently, 
in the Newtonian case a lot of matter takes up a relatively 
small amount of energy per particle, and the concept of DF 
is applicable, while in the case of MOND a small amount of 
matter takes up a lot of energy per particle, and DF is not 
an appropriate concept. The prolonged merging timescale in 
MOND is a consequence of the need to completely transform 
the orbits of the stars that are absorbing the original orbital 
energy. 



6 CONCLUSIONS 

Our fully non-linear simulations of the MOND dynamics of 
realistic stellar systems nicely confirm the analytical predic- 
tions of CB04. In particular, the simulations confirm the pre- 
dicted scaling of the DF timescale with the parameter 1 + TZ 
that measures the extent to which the actual acceleration 
exceeds that generated by the baryons alone in Newtonian 
gravity. Consequently, the astrophysical consequences listed 
in CB04 are confirmed. 

Prominent among these was the implication of short 
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DF timescales for the existence of globular clusters in dwarf 
spheroidal galaxies. Even simple Newtonian models of dwarf 
spheroidals and dwarf ellipticals predict that many observed 
globular clusters should spiral to the galaxy centr es within 
a Hubble time (|Tremainel [iQT^ : iLotz et al.|[200ll ). There- 
fore the existence of objects such as the Fornax dwarf 
spheroidal, which has five globular clusters but no stellar 
nucleus that could be the remains of clusters that have spi- 
ralled inwards, is a challenge to both Newtonian gravity and 
MOND. In the context of the standard cold-dark-matter 
theory, the survival of Fornax's globular cluster against 
DF is problematic because it can be explained only if the 
dark matter halo has an extended co re, and not a cen- 
tral cusp as predicted by the theory llGoerdt et ahl l2006l : 
ISanchez-Salcedo. Reves-Iturbide fc Hernandejl200^ From 
the point of view of MO ND, dwarf s pheroidals are in the 
deep-MOND regime fe.g. ICerhard fc Spcrgcl 1992), so the 
MOND inspiralling time of their globular clusters is as short 
as the dynamical time. Sanchez- Salcedo et al. (2006) eval- 
uated in detail the predictions of MOND and Newtonian 
gravity for the fate of the globular clusters of the Fornax 
dwarf spheroidal galaxy, and from the results of CB04 con- 
cluded that Fornax is much more problematic for MOND 
than for dark matter models. This conclusion gains strength 
from our confirmation of the central result of CB04. 
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